Introduction

A differential equation is an equation that involes the derivatives of a function. In this example, we make use of ordinary linear differential equations and discuss the analytic and numerical methods. They are in numerous fields like physics, biology, engineering and economics.

An Analytic solution is written below:

$$ \frac{dk}{dt} = ak$$ $$ \frac{dk}{k} = {a\cdot dt}$$ $$ \Rightarrow \int_{k_0}^{k} \frac{dk}{k} = \int_{0}^{t} a\cdot dt $$ $$ \Rightarrow ln | \frac{k}{k_0}| = a \cdot t$$ $$ \Rightarrow \frac{k}{k_0} = e^{a\cdot t} $$ $$ \Rightarrow k_0 \cdot e^{a\cdot t} $$ This shows exponential growth.

Another is written below:

$$ \frac{dI}{dt}= -d $$ $$ \Rightarrow \int_{I_0}^{I} {dI} = \int_{0}^{t} {-d}\cdot dt $$ $$ \Rightarrow I - I_0 = -d\cdot t $$ $$ I(t) = I_0 - dt $$

If $ c= a-b $ when a > b: $$\frac {dk}{dt} = c\cdot k$$ $$ k = k_0\cdot e^{-c\cdot t}$$ This showns exponential decay.

As differential equations become more complex, the tendency to use pen and paper to solve them will decrease.This is because of non-linearity where there will be no closed-form solutions and cannot be expressed as elementary functions (6). These equations will make use of iterative methods to find approximate solutions which would be more efficient to perform with code with the ability to act with other aspects of the equation much easily. It is also imperative to note that the solutions produced are approximate alongside the fact that analytical solutions can assume simple boundary conditions but numerical methods done on a computer can hander more realistic and complicated conditions allowing for us to change initial states. We are also able to perform graphical representations for better interpretation and analysis to understand trends.

Below is an example of how to solve a differential equation using the scipy library. We define a simple model: $$ \frac{dK}{dt} = a\cdot K - b\cdot K\cdot I $$ $$ \frac{dI}{dt} = c\cdot K\cdot I - d\cdot I $$ where K is capital, I is innovation. a is the growth rate of capital, b is the shrink rate of capital, c is the growth rate of innovation and d is the shrink rate of innovation.

This code block imports all the libaries with the last line of code used to suppress any errors when plotting the graphs after solving each differential equation iteratively.

import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from handcalcs.decorator import handcalc
import seaborn as sns
import logging # Suppress warnings 

logging.getLogger('scipy').setLevel(logging.ERROR)

This code block initalizes the parameters and defines the timestep and time series over which the code will run over.

a,b,c,d = 1,1,1,1
t0 = 0
tf = 30 #years
time_steps = int(tf / 0.01)
time_series = np.linspace(t0, tf, time_steps)

This code block translates the model written above in Latex into Python.

def my_model(t, y, a, b, c, d):
    K, I = y
    return [a*K - b*K*I, c*K*I - d*I]

This code block creates a helper function which takes in paramter K0 and I0 which are the initial states for Capital and Innovation. We make use of the solve_ivp to then solve the differential equation using the initial states.

def helper(K0, I0):
    sol = solve_ivp(fun = my_model, t_span = [t0, tf], y0 = [K0,I0], args=(a,b,c,d),t_eval = time_series)
    return sol

This code block now defines what I0 and K0 are. I0 is kept as a constant but K0 are a range of values with interval 2. The for loop is used to iterate through each value of K0 solving the equation and appending the solution to the array. Then, we select the second solution in the sols array to see if we can find any dynamics.

I0 = 0.1
K0 = [i for i in range(0,10,2)]
sols = []
for K0_ in K0:
    sols.append(helper(K0_, I0))

sol = sols[2]

We then plot a time series graph to see the variation between capital and innovation.

plt.plot(sol.t, sol.y[0], label = "Capital")
plt.plot(sol.t, sol.y[1], label = "Innovation")
plt.xlabel("Time")
plt.ylabel("Z")
plt.legend()

<matplotlib.legend.Legend at 0x7fd0d322dc40>

We then do the same but instead of a time series plot, we shall plot a phase diagram. A phase diagram is defined as a graphical representation of the trajectories x(t), y(t) vary against one another.

This code block is able to use the zip function to group the sols and K0 as one to then be able to plot the solutions against one another for both capital and innovation i.e the coupled pair of differential equations with the label being K_0. This is to plot it for each K0 as this is what is being varied.

for sol,K_0 in zip(sols,K0):
    plt.plot(sol.y[0], sol.y[1], label = K_0)
    plt.xlabel("Capital")
    plt.ylabel("Innovation")
plt.legend()

<matplotlib.legend.Legend at 0x7fd105d44770>

Now, we have not noticed anything interesting from this so we can introduce numerous other features and variables with the aspiration to obtain more dynamic observations.

$$ \frac{dK}{dt} = aK - bKIL + eLT - fRKL $$ $$ \frac{dI}{dt} = cKIL - dI + gTL $$ $$ \frac{dL}{dt} = -iL + jK $$ $$ \frac{dT}{dt} = kLKIR- oT + LfRK$$ $$ \frac{dR}{dt} = -mRLK + nLKI$$

The new terms introduced are L,T,R which are labour, technology and resources respectively.

For the first equation, we multiply capital by a constant as this is growth of capital and as capital grows, then capital itself shall grow with a key example being the stock market. Labour and technology are multipled by a constant as well because we are able to get capital as through the use of technology and labour(4), more things can be developed so greater capital. Labour is represented by people in this scenario. We then subtract by a constant mutlipled by K * I * L because Innovation is dependent on labour as people are needed for idea generation. And we multiply by capital because innovation is relative to capital in the sense that the more capital then the more innovation there shall be and therefore more capital will be used up. We then subtract by R * K * L multipled by a constant because Resources and Labour will use up more capital and this will again be proportional to the capital as we multiply it.

For the second equation, we multiply capital by innovation and labour because the more capital there is with greater labour then there is more people with greater money for idea generation to occur so innovation shall increase. The same will apply as technology is multiplied by labour as this causes idea generation (4)and therefore innovation is produced. At one point, old ideas will be deemed useless so we subtract a constant d mutliplied by innovation for this to occur.

For the third equation, capital is needed as the rate of capital can increase labour because more workers would mean that more money will be required(1). We then subtract this by a constant multipled by labour as these workers will retire soon and therefore the count will decrease.

For the fourth equation, Labour multiplied by capital, innovation and resources allow for technology to be made. More people with correct amount of money and ideas can create new and sustainable technology. This f * R * K is essentially converting the capital to technology through the utilisation of the resources. Then we subtract by a constant multiplied by T because technology can be deemed useless and old as new technology is developed.

For the fifth equation, Labour multiplied by Capital and innovation will allow for more resources to be produced as more people with greater money and innovation will drive the use of resources and will require more to be needed and we subtract it by L * K because more people with more capital relative to the amount of resources will cause the shrink of resources unfortunately.

Please see the table below for the input and outputs of the model:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.graph_objects as go
from tqdm.auto import tqdm
import pickle

This model is following concepts that have been written before but it now running upto 10000 times plotting a graph each time. We create a pickle file alongside a plot file to see how each parameter varies the plot to create interesting dynamics.

#This defines the model with the given set of equations above.
def my_model(t, y):
    K, I, L, T, R = y
    
    
    dK_dt = (a * K - b * K * I*L + e * L * T - f * R* K*L)
    dI_dt = (c * K * L * I - d * I + g * T * L) 
    dL_dt = (-i *L + j * K) 
    dT_dt = (k * L * K * I *R + L* f*K*R - o*T) 
    dR_dt = (-m * R * L * K + n *L*K*I) 
    return [dK_dt, dI_dt, dL_dt, dT_dt, dR_dt]
#The helper function solves the equations.
def helper(K0, I0):
    initial_conditions = [K0, I0, 0.1, 0.1, 0.1]  # Initial conditions for K, I, L, T, R
    sol = solve_ivp(fun=my_model, t_span=[t0, tf], y0=initial_conditions, t_eval=time_series, method = "LSODA")
    return sol
    

""""
for z in tqdm(range(100)):    
    t0 = 0
    tf = 100 # years
    time_steps = int(tf / 0.1)
    time_series = np.linspace(t0, tf, time_steps)
    
    #L, T, R are Labour, Technology and Resources respectively
    j = 0
    a = np.random.uniform(0,2) # We make use of the random function using uniform distribution to get a value between 0 and 2 for all the parameters. 
    b =  np.random.uniform(0,2) 
    c =  np.random.uniform(0,2)
    d = np.random.uniform(0,2)
    e = np.random.uniform(0,2)
    f =  np.random.uniform(0,2)
    g = np.random.uniform(0,2)
    h = np.random.uniform(0,2)
    i = np.random.uniform(0,2)
    j = np.random.uniform(0,2)
    k = np.random.uniform(0,2)
    l = np.random.uniform(0,2)
    m = np.random.uniform(0,2)
    n = np.random.uniform(0,2)
    o = np.random.uniform(0,2)
    arr = [a,b,d,e,f,g,h,i,j,k,l,m,n,o]

    
    I0 = 0.1
    K0_values = [i for i in range(2,3)]
    
    sols = []
    for K0 in K0_values:
        sols.append(helper(K0, I0))
#We now reorganise the solutions for Time, Capital, Innovation, Labour, Technology and resources into a dataframe such that it is now available as  a plot in a seaborn friendly format. 
    all_dfs = []
    for sol, K0 in zip(sols, K0_values):
        t_dataframe = pd.DataFrame(sol.t, columns=["Time"])
        K_dataframe = pd.DataFrame(sol.y[0], columns=["Capital"])
        I_dataframe = pd.DataFrame(sol.y[1], columns=["Innovation"])
        L_dataframe = pd.DataFrame(sol.y[2], columns=["Labor"])
        T_dataframe = pd.DataFrame(sol.y[3], columns=["Technology"])
        R_dataframe = pd.DataFrame(sol.y[4], columns=["Resources"])
        sol_df = pd.concat([t_dataframe, K_dataframe, I_dataframe, L_dataframe, T_dataframe, R_dataframe], axis=1) #1 is column
        sol_df["K0"] = K0
        all_dfs.append(sol_df)
    all_dfs = pd.concat(all_dfs)
    
    fig, axes = plt.subplots(3, 2, figsize=(20, 18))
    
    # Line plot for Capital
    sns.lineplot(x="Time", y="Capital", hue="K0", data=all_dfs, ax=axes[0, 0])
    axes[0, 0].set_title("Capital vs Time with varying K0 (Line Plot)")
    axes[0, 0].set_xlabel("Time")
    axes[0, 0].set_ylabel("Capital")
    axes[0, 0].grid(True)
    
    sns.lineplot(x="Time", y="Innovation", hue="K0", data=all_dfs, ax=axes[0, 1])
    axes[0, 1].set_title("Innovation vs Time with varying K0 (Line Plot)")    axes[0, 1].set_xlabel("Time")
    axes[0, 1].set_ylabel("Innovation")
    axes[0, 1].grid(True)
    
    sns.lineplot(x="Time", y="Labor", hue="K0", data=all_dfs, ax=axes[1, 0])
    axes[1, 0].set_title("Labor vs Time with varying K0 (Line Plot)")
    axes[1, 0].set_xlabel("Time")
    axes[1, 0].set_ylabel("Labor")
    axes[1, 0].grid(True)
    
    sns.lineplot(x="Time", y="Technology", hue="K0", data=all_dfs, ax=axes[1, 1])
    axes[1, 1].set_title("Technology vs Time with varying K0 (Line Plot)")
    axes[1, 1].set_xlabel("Time")
    axes[1, 1].set_ylabel("Technology")
    axes[1, 1].grid(True)
    
    sns.lineplot(x="Time", y="Resources", hue="K0", data=all_dfs, ax=axes[2, 0])
    axes[2, 0].set_title("Resources vs Time with varying K0 (Line Plot)")
    axes[2, 0].set_xlabel("Time")
    axes[2, 0].set_ylabel("Resources")
    axes[2, 0].grid(True)
    
    plt.tight_layout()
    plt.savefig(f"new_{z}.png")
#The pickle file has the parameters that cause the plot to occur. We store this file for parameter selection later.
    with open(f"new_{z}.pickle", 'wb') as handle: 
        pickle.dump(arr, handle, protocol=pickle.HIGHEST_PROTOCOL)
    plt.close()
"""

'"\nfor z in tqdm(range(100)):    \n    t0 = 0\n    tf = 100 # years\n    time_steps = int(tf / 0.1)\n    time_series = np.linspace(t0, tf, time_steps)\n    \n    #L, T, R are Labour, Technology and Resources respectively\n    j = 0\n    a = np.random.uniform(0,2) # We make use of the random function using uniform distribution to get a value between 0 and 2 for all the parameters. \n    b =  np.random.uniform(0,2) \n    c =  np.random.uniform(0,2)\n    d = np.random.uniform(0,2)\n    e = np.random.uniform(0,2)\n    f =  np.random.uniform(0,2)\n    g = np.random.uniform(0,2)\n    h = np.random.uniform(0,2)\n    i = np.random.uniform(0,2)\n    j = np.random.uniform(0,2)\n    k = np.random.uniform(0,2)\n    l = np.random.uniform(0,2)\n    m = np.random.uniform(0,2)\n    n = np.random.uniform(0,2)\n    o = np.random.uniform(0,2)\n    arr = [a,b,d,e,f,g,h,i,j,k,l,m,n,o]\n\n    \n    I0 = 0.1\n    K0_values = [i for i in range(2,3)]\n    \n    sols = []\n    for K0 in K0_values:\n        sols.append(helper(K0, I0))\n#We now reorganise the solutions for Time, Capital, Innovation, Labour, Technology and resources into a dataframe such that it is now available as  a plot in a seaborn friendly format. \n    all_dfs = []\n    for sol, K0 in zip(sols, K0_values):\n        t_dataframe = pd.DataFrame(sol.t, columns=["Time"])\n        K_dataframe = pd.DataFrame(sol.y[0], columns=["Capital"])\n        I_dataframe = pd.DataFrame(sol.y[1], columns=["Innovation"])\n        L_dataframe = pd.DataFrame(sol.y[2], columns=["Labor"])\n        T_dataframe = pd.DataFrame(sol.y[3], columns=["Technology"])\n        R_dataframe = pd.DataFrame(sol.y[4], columns=["Resources"])\n        sol_df = pd.concat([t_dataframe, K_dataframe, I_dataframe, L_dataframe, T_dataframe, R_dataframe], axis=1) #1 is column\n        sol_df["K0"] = K0\n        all_dfs.append(sol_df)\n    all_dfs = pd.concat(all_dfs)\n    \n    fig, axes = plt.subplots(3, 2, figsize=(20, 18))\n    \n    # Line plot for Capital\n    sns.lineplot(x="Time", y="Capital", hue="K0", data=all_dfs, ax=axes[0, 0])\n    axes[0, 0].set_title("Capital vs Time with varying K0 (Line Plot)")\n    axes[0, 0].set_xlabel("Time")\n    axes[0, 0].set_ylabel("Capital")\n    axes[0, 0].grid(True)\n    \n    sns.lineplot(x="Time", y="Innovation", hue="K0", data=all_dfs, ax=axes[0, 1])\n    axes[0, 1].set_title("Innovation vs Time with varying K0 (Line Plot)")    axes[0, 1].set_xlabel("Time")\n    axes[0, 1].set_ylabel("Innovation")\n    axes[0, 1].grid(True)\n    \n    sns.lineplot(x="Time", y="Labor", hue="K0", data=all_dfs, ax=axes[1, 0])\n    axes[1, 0].set_title("Labor vs Time with varying K0 (Line Plot)")\n    axes[1, 0].set_xlabel("Time")\n    axes[1, 0].set_ylabel("Labor")\n    axes[1, 0].grid(True)\n    \n    sns.lineplot(x="Time", y="Technology", hue="K0", data=all_dfs, ax=axes[1, 1])\n    axes[1, 1].set_title("Technology vs Time with varying K0 (Line Plot)")\n    axes[1, 1].set_xlabel("Time")\n    axes[1, 1].set_ylabel("Technology")\n    axes[1, 1].grid(True)\n    \n    sns.lineplot(x="Time", y="Resources", hue="K0", data=all_dfs, ax=axes[2, 0])\n    axes[2, 0].set_title("Resources vs Time with varying K0 (Line Plot)")\n    axes[2, 0].set_xlabel("Time")\n    axes[2, 0].set_ylabel("Resources")\n    axes[2, 0].grid(True)\n    \n    plt.tight_layout()\n    plt.savefig(f"new_{z}.png")\n#The pickle file has the parameters that cause the plot to occur. We store this file for parameter selection later.\n    with open(f"new_{z}.pickle", \'wb\') as handle: \n        pickle.dump(arr, handle, protocol=pickle.HIGHEST_PROTOCOL)\n    plt.close()\n'

Above, we have analysed the images and we have seen that the plot "new_13.pickle" has a relaxation oscillation with a dynamic that can be explored in further detail compared to the other plots which were less dynamic. This is the plot "new_13.png" for visulisation purposes:

This section of code is run 14 times i.e the number of parameters and what is being done is that we multiply each of the parammters by two separately to see what plot is being produced.

for z_ in range(14):
    
    with open('new_13.pickle', 'rb') as handle: 
        arr_loaded = pickle.load(handle)
        arr_loaded[z_] = arr_loaded[z_] * 10 # Multiplying the parameters one at a time by 10
    
    t0 = 0
    tf = 100  # years
    time_steps = int(tf / 0.1)
    time_series = np.linspace(t0, tf, time_steps)
    
    
    a,b,d,e,f,g,h,i,j,k,l,m,n,o = arr_loaded
    
    
    I0 = 0.1
    K0_values = [i for i in range(2,3)]
    
    sols = []
    for K0 in K0_values:
        sols.append(helper(K0, I0))
    
    all_dfs = []
    for sol, K0 in zip(sols, K0_values):
        t_dataframe = pd.DataFrame(sol.t, columns=["Time"])
        K_dataframe = pd.DataFrame(sol.y[0], columns=["Capital"])
        I_dataframe = pd.DataFrame(sol.y[1], columns=["Innovation"])
        L_dataframe = pd.DataFrame(sol.y[2], columns=["Labor"])
        T_dataframe = pd.DataFrame(sol.y[3], columns=["Technology"])
        R_dataframe = pd.DataFrame(sol.y[4], columns=["Resources"])
        sol_df = pd.concat([t_dataframe, K_dataframe, I_dataframe, L_dataframe, T_dataframe, R_dataframe], axis=1) #1 is column
        sol_df["K0"] = K0
        all_dfs.append(sol_df)
        all_dfs = pd.concat(all_dfs)
        
        fig, axes = plt.subplots(3, 2, figsize=(20, 18))
    
    # Line plot for Capital
    sns.lineplot(x="Time", y="Capital", hue="K0", data=all_dfs, ax=axes[0, 0])
    axes[0, 0].set_title("Capital vs Time with varying K0 (Line Plot)")
    axes[0, 0].set_xlabel("Time")
    axes[0, 0].set_ylabel("Capital")
    axes[0, 0].grid(True)
    
    sns.lineplot(x="Time", y="Innovation", hue="K0", data=all_dfs, ax=axes[0, 1])
    axes[0, 1].set_title("Innovation vs Time with varying K0 (Line Plot)")
    axes[0, 1].set_xlabel("Time")
    axes[0, 1].set_ylabel("Innovation")
    axes[0, 1].grid(True)
    
    sns.lineplot(x="Time", y="Labor", hue="K0", data=all_dfs, ax=axes[1, 0])
    axes[1, 0].set_title("Labor vs Time with varying K0 (Line Plot)")
    axes[1, 0].set_xlabel("Time")
    axes[1, 0].set_ylabel("Labor")
    axes[1, 0].grid(True)
    
    sns.lineplot(x="Time", y="Technology", hue="K0", data=all_dfs, ax=axes[1, 1])
    axes[1, 1].set_title("Technology vs Time with varying K0 (Line Plot)")
    axes[1, 1].set_xlabel("Time")
    axes[1, 1].set_ylabel("Technology")
    axes[1, 1].grid(True)
    
    sns.lineplot(x="Time", y="Resources", hue="K0", data=all_dfs, ax=axes[2, 0])
    axes[2, 0].set_title("Resources vs Time with varying K0 (Line Plot)")
    axes[2, 0].set_xlabel("Time")
    axes[2, 0].set_ylabel("Resources")
    axes[2, 0].grid(True)
    plt.suptitle(z_) # Which parameter number is being varied.
    plt.tight_layout()
    plt.show()
    fig = plt.figure(figsize=(12, 10))
    ax = fig.add_subplot(111, projection='3d')
    for K0 in K0_values:
        df = all_dfs[all_dfs["K0"] == K0]
        sc = ax.scatter(df["Capital"], df["Innovation"], df["Labor"], c=df["Technology"], cmap='viridis', s=df["Resources"] * 10, label=f'K0={K0}')
    ax.set_xlabel('Capital')
    ax.set_ylabel('Innovation')
    ax.set_zlabel('Labor')
    ax.set_title('Phase Plot of Capital, Innovation, Labor with Technology (color) and Resources (size)')
    cbar = fig.colorbar(sc, ax=ax)
    cbar.set_label('Technology')
    ax.legend()
    plt.show()

    sns.set(style="darkgrid")

    # Define pairs of variables to plot
    variable_pairs = [('Capital', 'Labor'), ('Capital', 'Innovation'), ('Capital', 'Technology'),
                      ('Labor', 'Innovation'), ('Labor', 'Technology'), ('Innovation', 'Technology')]
    
    # Plotting each pair
    plt.figure(figsize=(15, 10))
    for i, (var1, var2) in enumerate(variable_pairs, 1):
        plt.subplot(3, 2, i)
        sns.scatterplot(x=var1, y=var2, hue="K0", data=all_dfs)
        plt.title(f"{var1} vs {var2}")
        plt.xlabel(var1)
        plt.ylabel(var2)
        plt.legend(title='K0', loc='upper right')
    
    plt.tight_layout()
    plt.show()
        

 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.1851401400309D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.1851401400309D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.1851401400309D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.1586960318103D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.1586960318103D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.1586960318103D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.1231698297079D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.1231698297079D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.1231698297079D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.2844551874787D+01   r2 =  0.9779798092547D-16
 lsoda--  above warning has been issued i1 times.  
       it will not be issued again for this problem
      in above message,  i1 =        10
 lsoda--  at t (=r1), too much accuracy requested  
       for precision of machine..  see tolsf (=r2) 
      in above,  r1 =  0.2844551874787D+01   r2 =                  NaN

/tmp/ipykernel_270099/1837334885.py:9: RuntimeWarning: overflow encountered in scalar multiply
  dT_dt = (k * L * K * I *R + L* f*K*R - o*T)
/tmp/ipykernel_270099/1837334885.py:9: RuntimeWarning: invalid value encountered in scalar subtract
  dT_dt = (k * L * K * I *R + L* f*K*R - o*T)
/tmp/ipykernel_270099/1837334885.py:6: RuntimeWarning: invalid value encountered in scalar subtract
  dK_dt = (a * K - b * K * I*L + e * L * T - f * R* K*L)
/tmp/ipykernel_270099/1837334885.py:7: RuntimeWarning: invalid value encountered in scalar subtract
  dI_dt = (c * K * L * I - d * I + g * T * L)
/tmp/ipykernel_270099/1837334885.py:10: RuntimeWarning: invalid value encountered in scalar add
  dR_dt = (-m * R * L * K + n *L*K*I)
/home/sanjhay/miniconda3/envs/test-evironment/lib/python3.12/site-packages/scipy/integrate/_ivp/lsoda.py:161: UserWarning: lsoda: Excess accuracy requested (tolerances too small).
  solver._y, solver.t = integrator.run(

 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.3479128743293D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.3479128743293D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.2933765687164D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.2933765687164D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.2933765687164D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.2662430602652D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.2662430602652D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.2662430602652D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.1977675397648D-15
 lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.6380837225686D+01   r2 =  0.1977675397648D-15
 lsoda--  above warning has been issued i1 times.  
       it will not be issued again for this problem
      in above message,  i1 =        10
 lsoda--  at t (=r1), too much accuracy requested  
       for precision of machine..  see tolsf (=r2) 
      in above,  r1 =  0.6380837225686D+01   r2 =                  NaN

/tmp/ipykernel_270099/1837334885.py:9: RuntimeWarning: overflow encountered in scalar multiply
  dT_dt = (k * L * K * I *R + L* f*K*R - o*T)
/tmp/ipykernel_270099/1837334885.py:9: RuntimeWarning: invalid value encountered in scalar subtract
  dT_dt = (k * L * K * I *R + L* f*K*R - o*T)
/tmp/ipykernel_270099/1837334885.py:6: RuntimeWarning: invalid value encountered in scalar subtract
  dK_dt = (a * K - b * K * I*L + e * L * T - f * R* K*L)
/tmp/ipykernel_270099/1837334885.py:7: RuntimeWarning: invalid value encountered in scalar subtract
  dI_dt = (c * K * L * I - d * I + g * T * L)
/tmp/ipykernel_270099/1837334885.py:10: RuntimeWarning: invalid value encountered in scalar add
  dR_dt = (-m * R * L * K + n *L*K*I)
/home/sanjhay/miniconda3/envs/test-evironment/lib/python3.12/site-packages/scipy/integrate/_ivp/lsoda.py:161: UserWarning: lsoda: Excess accuracy requested (tolerances too small).
  solver._y, solver.t = integrator.run(

Above, there are three sets of plots for varying each constant. There are some interesting features to note like changing the constants in the labour equation, technology equation and innovation equation. Now, I have decided to pick out important plots of the model to attempt to analyse them.

A general trend seen in the 2D Phase plots for the variable pairs are that for Capital vs Technology, there is a suggestion of oscillatory behaviour and that the system can undergo cycles where technology increases with capital, reaches a peak, and then declines as capital continues to rise or fall. The closed loop patterns suggest stable limit cycles where ths sytem returns to the same state. The Labour vs Innovation indicates a more complex dynamic where labour and innovation affect each other in a non-linear manner suggesting bistable states as a decrease in innovation causes an increase in labour. The Labour vs Technology represents regular cycles of increasing and decreasing technology as labour varies with a more emphasised amplitude as opposed to other plots seen. Lastly, the Innovation vs Technology has distinct oscillations that are less periodic but also more complex and shows how innovation drives technology, vice versa.

At the lower values of Capital and Labour, the system starts with darker colours with resources as relatively small. As the system evolves, both labour and capital increase and there is an outcome of greater technology indicated by the lighter colours showing how technology is increasing rapidly whilst innovation increases. Soon, the system near the end stabilizes into repetitive cycles with periodic variations for technology, capital, labour and innovation.

We can see a general trend in this graph where there are dynamics to be noted. There are successive periodic oscillations which are constantly generated for Capital, Innovation. Howeveer, for Labour it starts from 0 and then as time increases, the graph never goes back down to 0 but decreases slightly which is interesting to note. For technology, it had a boom before suddenly crashing down and then increasing. It is important to note that it increases and then decreases producing a passive oscillation.

Conclusion

It is imperative to note the effect seen as the parameters are multiplied by 10 one by one keeping the rest constant. It gives us an insight of how changing one parameter can change the entire model. In the future, we can try and incorporate more stochasticity i.e noise whilst also increasing the complexity of the model, perchance using partial derivatives. This model has outlined the impacts that capital, technology, resources, innovation and labour have against one another, how they progress over time and lastly how they interact with one another.

Bibleograhy

1. https://futuresofwork.co.uk/2022/11/29/technologys-impact-on-the-labour-market-different-this-time/
2. https://core.ac.uk/download/pdf/82214475.pdf
3. https://www.oreilly.com/library/view/numerical-methods/9789332503465/xhtml/chapter009.xhtml4.
4. https://rcc.harvard.edu/knowledge-technology-and-complexity-economic-growth#:~:text=Technological%20progress%20allows%20for%20the,used%20in%20production%20are%20complex.
5. https://en.wikipedia.org/wiki/Euler_method
6. https://x-engineer.org/euler-integration/

I have referred to this using in text-citations and the corresponding numbers.